/*!  \file  SteVD.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      07/07/2013 12:01:05 PM
 *   \copyright GNU Public License.
 */
#include "SteVD.h"
#include "Matrix.h"
#include "setThreads.h"
extern "C" {
void  dstevd_(
const char * jobz, 
const int  *    n,
double *  d, 
double *  e, 
double *  z, 
const int *   ldz, 
double * work, 
const int * lwork, 
int * iwork,
const int * liwork, 
int * info );
}
SteVD::SteVD(): TriEigen("Lapack dstevd") {	
	setCompileTime(__DATE__, __TIME__); 
	canBeParallel = true;
}
bool SteVD::setNumberOfThreads(int nThreads_) {	
	assert(canBeParallel);
	SET_NUM_THREADS(nThreads_);
	nThreads = nThreads;
	return true;
}

void SteVD::compute(const Matrix & T, Matrix & Sigma, Matrix & Q ){
	reset();
	assert(T.isTriDiagonal() && Sigma.isDiagonal() && Q.isOrthogonal());
	char jobz = 'V';	int n = T.cols(), info = -1, ldq = Q.lda();
	double * e = NULL; allocate(e, n); 
	double * d    = Sigma.getModData();
	double * q    = Q.getModData(); Q.makeColMajor();
	memcpy(d, T.getData()    ,    n * sizeof(*d));
	memcpy(e, T.getOffDiagonal(), (n-1)* sizeof(*d));
	int lwork = -1, liwork = -1;
	double optW; int optIW;
	// calculates the optimal sizes of the WORK and IWORK arrays
	dstevd_(&jobz, &n, d, e, q,&ldq, &optW, &lwork, &optIW, &liwork, &info);
	success &= info == 0;
	if (info == 0) {
		lwork = (int) optW; 	liwork = (int) optIW;
		double *work = NULL; int *iwork = NULL;
		allocate(work, lwork); allocate(iwork, liwork);
		if (work && iwork ) {
			dstevd_(&jobz,&n, d, e, q, &ldq, work, &lwork, iwork, &liwork, &info);
			success &= info == 0;
		}
		delocate (iwork);	delocate ( work);
	}
	delocate(e);
};
BiNESteVD::BiNESteVD() : BiSVD( "NE dstevd" ) {
	aUpper = true; 
	aLower = false; 
	aUColMajor = true; 
	aVColMajor = true; 
	VColVec = false;
	setTriEigen(*(new SteVD()));
	setCompileTime(__DATE__, __TIME__); 
};
BiNE2SteVD::BiNE2SteVD() : BiSVD( "NE2 dstevd" ) {
	aUpper = true; 
	aLower = false; 
	aUColMajor = true; 
	aVColMajor = true; 
	VColVec = false;
	setTriEigen(*(new SteVD()));
	setCompileTime(__DATE__, __TIME__);
};
//REGISTER(SteVD);
//REGISTER(BiNESteVD);
//REGISTER(BiNE2SteVD);
